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ABSTRACT 

Zonal and meridional axisymmetric flows can deeply impact the rotational and chem¬ 
ical evolution of stars. Therefore, momentum exchanges between waves propagating 
in stars, differential rotation, and meridional circulation must be carefully evaluated. 
In this work, we study axisymmetric mean flows in rapidly and initially uniformly 
rotating massive stars driven by small amplitude non-axisymmetric tt-driven oscilla¬ 
tions. We treat them as perturbations of second-order of the oscillation amplitudes 
and derive their governing equations as a set of coupled linear ordinary differential 
equations. This allows us to compute 2-D zonal and meridional mean flows driven 
by low frequency g- and r-modes in slowly pulsating B (SPB) stars and p-modes in 
j3 Cephei stars. Oscillation-driven mean flows usually have large amplitudes only in 
the surface layers. In addition, the kinetic energy of the induced 2-D zonal rotational 
motions is much larger than that of the meridional motions. In some cases, meridional 
flows have a complex radial and latitudinal structure. We find pulsation-driven and 
rotation-driven meridional flows can have similar amplitudes. These results show the 
importance of taking wave - mean flow interactions into account when studying the 
evolution of massive stars. 

Key words: hydrodynamics - waves - stars: rotation - stars: oscillations - stars: 
evolution - stars: massive 


1 INTRODUCTION 


It is well known that differential rotation and related viscous turbulent transport, structural adjustments, and applied torques 
at stellar surfaces drive large-scale meridional circulation in stellar radiation zones (e.g., Zahn 1992; Rieutord 2006; Decressin 
et al. 2009; Hypolite & Rieutord 2014). When magnetic fields and internal gravity waves are neglected, the magnitudes of the 
radial component of such meridional flow ( v r ) scales as (e.g., Kippenhahn et al. 2012) 


Vr:MC 


LR 2 Sr 
GM 2 2t rGp’ 


(1) 


where L, M, R, p and Q are the luminosity, mass, radius, mean density and angular velocity of the star respectively. It is 
anticipated that such meridional circulation in rotating stars mix material and transport angular momentum in their interior 
during their evolution (e.g., Zahn 1992; Talon et al. 1997; Maeder & Zahn 1998; Meynet & Maeder 2000; Mathis & Zahn 
2004). 

Stellar non-radial pulsations in rotating stars can also drive such non-oscillatory, large-scale fluid motions, which we 
may call mean flows. Because of dissipative mechanisms, co-rotation resonances and breaking mechanisms (e.g., Andrews & 
McIntyre 1978a; Lindzen 1981; Goldreich & Nicholson 1989), pulsations transport angular momentum, which induces mean 
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zonal flows leading to differential rotation. Indeed, angular momentum transport by internal gravity waves takes place in 
rotating stars, if the wave propagation is accompanied by thermal diffusion/critical layers/non-linear breaking processes that 
lead to damping and/or excitation of the waves (e.g., Press 1981; Schatzman 1993, 1996; Zahn et al. 1997; Alvan et al. 
2013; Rogers et al. 2013). For example, if prograde internal gravity waves propagating in a rotating star damp at a certain 
place, they will deposit their angular momentum, which accelerates zonal rotational motion of the fluid there. On the other 
hand, when prograde waves are excited in a place, they extract angular momentum from the flow at this location, leading 
to deceleration of the rotational fluid motion. The opposite results are obtained for retrograde waves. Therefore, extraction 
and deposition of angular momentum from rotational fluid motion by waves change the zonal velocity field in the interior of 
rotating stars. Because of the related modification of angular velocity gradients, of the radial and latitudinal components of 
the momentum equation, and of the equation of heat, waves modify the general 3-D dynamical balance in stellar radiation 
zones. As a consequence, they also profoundly affect large-scale meridional flows in the radial and latitudinal directions in 
these regions (e.g., Bretherton 1969; Andrews & McIntyre 1978a; Holton 1982; Mathis et al. 2013; Belkacem et al. 2015a). 
In the case of small-amplitude linear pulsations, magnitudes of mean flows are of second-order of the amplitudes and their 
characteristic timescale is the growth or damping timescales of the waves. Then, these pulsation-driven mean flows will impact 
material mixing and angular momentum transport processes in rotating stars and will affect pulsations in return. 

These processes have been relatively well studied for the case of low-mass stars (e.g., Mathis et al. 2013), but poorly 
investigated for massive stars. Therefore, it is important to investigate the velocity fields and the magnitudes of both zonal 
and meridional mean flows driven by stellar pulsations in rotating stars and to investigate their geometrical properties and 
amplitude in the case of massive stars. 

The problem of angular momentum transport and redistribution by internal gravity waves in rotating stars has been 
investigated by many authors for cool stars. For example, for an explanation of the uniform rotation observationally inferred 
for the radiative core of the Sun until 0.2 Rq thanks to helioseismology (e.g., Schou et al. 1998; Garcia et al. 2007), Schatzman 
(1993) suggested that angular momentum transport by internal gravity waves generated in the external convective envelope 
can play an essential role to extract angular momentum [see also Zahn et al. (1997); Talon et al. (2002); Talon & Charbonnel 
(2005); Rogers et al. (2008); Mathis et al. (2008) in the case where the Coriolis acceleration is taken into account; Kumar 
et al. (1999); Rogers & MacGregor (2010); Mathis & de Brye (2012) in the magnetized case]. More recently, thanks to new 
asteroseismic information obtained on the weak differential rotation between the surface and the core of low-mass subgiant 
and red giant stars (e.g., Beck et al. 2012; Mosser et al. 2012; Deheuvels et al. 2012, 2014), the potential extraction of angular 
momentum by internal gravity waves and normal mixed gravito-acoustic modes have been studied (Talon & Charbonnel 2008; 
Fuller et al. 2014; Belkacem et al. 2015a,b). 

In the case of massive stars, the transport of angular momentum by internal gravity waves has been less studied (Lee & 
Saio 1993; Pantillon et al. 2007; Lee et al. 2014). Zahn (1975, 1977) investigated the case of gravity waves tidally excited in 
massive main-sequence stars by the orbital motion of a companion star in a binary system (see also Goldreich & Nicholson 
1989). It allowed him to estimate the timescales of circularization of the orbit and of synchronization of the stellar rotation 
with the orbital motion. However, the impact of internal gravity waves on large-scale meridional flows in these stars has never 
been studied. 

In this context, one important issue is the case of Be stars. Be stars are rapidly rotating main-sequence late-O, B, and 
early-A stars that have circumstellar gaseous discs, which are responsible for the generation of emission lines (see Rivinius et 
al. 2013 for a complete review). Gaseous discs in Be stars are thought to be viscous Keplerian discs (Lee et al. 1991). Although 
mechanisms for disc formation have not been completely identified yet, rapid rotation and pulsations of these stars must be 
key ingredients. Moreover, to sustain the discs for several decades, a good amount of angular momentum must be continuously 
supplied to them. Ando (1983, 1986) employed a wave mean flow interaction theory (e.g., Andrews & McIntyre 1976, 1978a,b; 
Dunkerton 1980; Grimshaw 1984; Craik 1988; Pedlosky 1982; Biihler 2014) to explain episodic mass loss phenomena observed 
in Be stars (see also Lee et al. 1991, 2014). As suggested by Lee (2013), if enough angular momentum is supplied to the surface 
regions of rapidly rotating stars, we can construct a steady system composed of a rotating star and a viscous Keplerian disc 
around the star. As a mechanism to provide the surface layers with angular momentum, we may invoke the rot at ion-driven 
meridional circulation (Meynet & Maeder 1997; Ekstrom et al. 2008; Granada et al. 2013) and/or internal gravity waves, 
which may be destabilized by the opacity mechanism or stochastically excited by convective motion in the core (see e.g., 
Rogers et al. 2013; Lee et al. 2014). 

In this paper, we are thus interested in velocity fields of pulsation driven mean zonal and meridional flows in rotating 
massive stars, and we investigate whether the velocity fields of the mean flows can be favorable to a disc formation around 
Be stars. Therefore, we calculate pulsation-driven mean flows in the interior of rotating SPB and /3 Cephei stars, in which 
pulsations are mainly excited by the iron opacity bump (k) mechanism. Assuming that the pulsation-driven velocity fields 
are of second-order of the linear pulsation amplitudes, we derive for the second-order perturbations a set of linear ordinary 
differential equations, which have inhomogeneous terms generated by non-linear products of the eigenfunctions of the pulsation 
mode. To derive analytical expressions of these non-linear terms, we introduce and use the so-called Spin Weighted Spherical 
Harmonics and their complex and tedious derivation is detailed in a complete appendix. The set of linear differential equations 
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we solve is given in §2. In §3, we calculate and discuss the properties of second-order velocity fields driven by unstable g- and 
r-modes in SPB stars and p-modes in f3 Cephei stars assuming moderate rotation of the stars. In §4, we summarize the results 
we obtained, and we discuss drawbacks and perspectives of this work. 


2 THEORETICAL FORMALISM FOR MEAN FLOWS INDUCED BY PULSATIONS 


2.1 Dynamical equations 


The basic equations governing fluid motions in a star are given by 


dv 1 

——b v ■ Vd = —Vp — VH?, 
at p 

(2) 

V 2 $ = 4nGp, 

( 3 ) 

l + v-<H = o, 

( 4 ) 


( 5 ) 

F = —AVT with A = 16 <tsbT 3 /3kp, 

(6) 


where v is the velocity vector of the fluid, p is the pressure, p is the mass density, T is the temperature, s is the specific 
entropy, F is the energy flux vector, <f> is the gravitational potential, e is the nuclear energy generation rate per gram, k is 
the opacity, G is the gravitational constant, and ctsb is the Stefan-Boltzmann constant. 

We assume that the star pulsates in small amplitudes around the hydrostatic and radiative equilibrium state. If pulsation 
amplitudes are sufficiently small, any physical quantity / (x, t) of the star may be represented by 

= / (0) (x) + f'(x, t) + / (2) (a;, t) + • • •, (7) 

where /® denotes equilibrium quantities, /' Eulerian perturbations of first-order, and /® Eulerian perturbations of second- 
order of the pulsation amplitude. The velocity field v (x, t) may also be given by 

v(x,t) = v ^ (x) + v'(x, t) + v^ (x, t) + - - - , (8) 

and the equilibrium state is assumed to be that of a uniformly rotating star so that, in spherical polar coordinates (r, #,</>), 

= rsin#fle^>, (9) 

where fl is the angular velocity of rotation and assumed to be constant, and e</, is the unit vector in the azimuthal direction. 
In this paper, we ignore rotational deformation of the equilibrium state so that depends only on the radial distance r from 
the center of the star. We also employ the Cowling approximation, neglecting the Euler perturbation of the gravitational 
potential <f>. 


2.2 First-order pulsation quantities 

For a given equilibrium state of a star, we solve the linear oscillation equation to obtain wave quantities represented by f'. 
For uniformly rotating stars, the oscillation equations, which we solve to obtain normal modes, are found, for example, in Lee 
& Saio (1987) and Lee & Baraffe (1995). Since separation of variables between the radial (r) and angular ( 8,(j >) coordinates 
is not possible for perturbations in rotating stars, we use finite series expansion to represent the perturbations in terms of 
spherical harmonic functions Y; m (#,(/>) (e.g., Lee & Saio 1987). For the Lagrangian displacement vector £, we may write 

£(*, t) = e 1CT *£(*) = e 1CTt (£ r e r + £ e e 0 + ^e^,), (10) 

where e r , eg, and e^ are the orthonormal basis vectors in spherical polar coordinates, a denotes the oscillation frequency 
observed in an inertial frame. For a given azimuthal order m, the components of £(cc) are given by 

Jmax 

Zr(x) = r'£S l flr)Y l ™(e,<l>), (11) 

l=i 

Jmax 

&(x) = 

3 = 1 

&(*): = r J2 

3 = 1 




( 13 ) 




( 12 ) 
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and the Eulerian pressure perturbation, p'(x,t) = p'(x)e lat , is given by 

3 max 

p>)=£p«i( r ) y «rM’ (i4) 

3=1 

where lj = 2 (j — 1 ) + \m\ and lj = lj + 1 for even modes, and lj = 2j — 1 + |m| and = lj — 1 for odd modes for 
= 1 ,2,3*--, jmax- The angular dependence of p'(x) is symmetric (antisymmetric) about the equator for even (odd) modes. 
In addition, we have 

v = 5v — £ ■ Vr® = -y- — £ ■ Vt/ 0 ^ = ka£, (15) 

at 

where ui = a + mfl is the oscillation frequency observed in the corotating frame of the star, and Sv denotes the Lagrangian 
perturbation of the velocity vector. 


2.3 Pulsation driven mean flows 

2.3.1 Expansion of second-order equations 


We use a theory of wave-mean flow interaction to discuss axisymmetric flows driven by non-axisymmetric oscillations in 
rotating stars. We regard the axisymmetric flows as mean flows, which contain both zero th and second-order contributions in 
the oscillation amplitude. The zero th order quantities /® are those of equilibrium state and are independent of time t. The 
second-order quantities /^ carry the time dependence of the mean flow. To derive governing equations for the second-order 
quantities f^ 2 \ we use the zonal averaging defined by 

(16) 

and if we assume 


/' = 0, (17) 

we have 


/ = /(°)+/( 2 ). 


(18) 


Here, we have ignored higher order terms with k ^ 3. Hereafter, we simply write /® and /^ for /(°) and /( 2 ). Because 
of the zonal averaging, / < '°- ) and /® are independent of tj>. 

When a non-axisymmetric oscillation mode is excited to attain a small but finite amplitude, non-oscillatory fluid flows 
may arise as a result of second-order effects of the oscillation. Applying the zonal averaging to the basic equations, we obtain 
a set of differential equations that govern second-order quantities: 


dv (2) 

dt 


+ v 


( 0 ) 


• V« (2) + v (2) • Vv (0) + 


J—\7p( 2 ) _l Q P _ - e 

9 (0)Vp +ff p(0 ) e r 


—v' ■ Vv' + 





^ + V • (p ( 0 V 2) ) = -V ■ {p'v'), 

F (2) = -A (0) VT (2) - A (2) VT (0) - A'VT', 


(19) 

( 20 ) 
( 21 ) 


p (0) T (0) + w (2) • Vs (0) ) = -p (0) T (0) + v' ■ Vs(°) + wC0) . V S ') - p (0) T (0 V ■ Vs' 

+p (0) e (2) + p (2) e (0) - V ■ F (2) + 7?, (22) 

where p = —(p^)~ 1 dp^ /dr = GM r /r 2 with M r = J//4nr' 2 p^ dr'. Because of the zonal averaging, we have for non- 
axisymmetric oscillations 


fg' = £ K (f(r, 0)e im ^) K (g\r, 0)e ira *) d(/> — (/'* p') = ^ (/'<?'*) , 

where the asterisk (*) indicates complex conjugation. 


(23) 


2.3.2 Solving the system of equations 

Since the second-order quantities are assumed axisymmetric, we expand the velocity perturbation i/ 2 - 1 using spherical harmonic 
functions Y t ° as 
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vfl\x,t) = Y v ( s} k {r,t) Y ° k {9,4>), 
k =1 

^max r\ 

= Y ^M) ™^°(M), 

k=1 

^max o 

v£\x,t) = ~Y d T?i'fl r ’ t )-Q(j Y ° k ( e ' ( t > ')' 


(24) 

(25) 

(26) 


and scalar physical quantities, such as pressure perturbation p' T> for example, as 

fcmax 

P (2) {x,t ) = Y Ph ( r >t)*l°(M). 

fc=l 


(27) 


where 4 = 2(fc — 1) and l'k =4 + 1 for fc = 1, 2, • • •, fc max . Note that for both even and odd linear modes the r and 
4> components of the right-hand-side of equation for example, are symmetric about the equator of the star and the 0 
component is antisymmetric. Since the expansion coefficients of v ® contain contributions from products Y / ^ rn Y" 1 , we have to 
set the expansion length fc max ~ 2 x ji max for the second-order coefficients. The time dependence is included in the expansion 
coefficient. 

By substituting the expansions (Hi to 03) into equations m to multiplying by a given spherical harmonic 

function, and integrating over spherical surface, we derive a finite set of differential equations for the expansion coefficients, 
which depends on r and t. To compute the integrals corresponding to non-linear terms, such as J(Y° )* (v' • Vo') r sin ddddcj), we 
have to evaluate angular integration of products of three spherical harmonic functions, which can be systematically carried out 
by introducing spin-weighted spherical harmonic functions sY™(6, (j>) , the definition and properties of which are summarized 
in Appendix A (see also, e.g., Newman & Penrose 1966; Varshalovich et al. 1988). In Appendix Bl, we introduce a new 
set of basis vectors e r , e q and e q , and rewrite the velocity and displacement vectors, and basic equations (1191) to (12211 on 
this basis. Then, we derive explicit expressions of non-linear terms by using spin-weighted spherical harmonics (Appendix 
B2). This allows us to obtain ordinary differential equations for the radial functions of the velocity field and scalar quantities 
(Appendixes B3 and B4). 

Using vector notation for the dependent variables of second-order defined as 



( vfl/rao\ 


( p[f/9 r P {0) \ 


/ L { r}jLr ) \ 

Zi = 

rfjflrao 

, *2 = 

o' 

, *3 = 

T (2) / r (0) 

L/ r,l 2 t- Lr 


l : ) 


1 : ) 


l / 





(y 



(pf\ 

Zh = 

V^hti (2) l2 /rao 

, z t = 


/■ rao 

(2) 

1 p = 

Pl2 


l : ) 


\ 

■ ) 


l : ) 


we rewrite equations (IB68I) to (IB71I) as 


dz 


I - A 1 / 2 


dr 

dZ2 


Zh — ( 3 + rA — ¥- 
1 i 


( 2 ) 


_ 9 P' 

1 dr p(°) 


H° 


crop 


( 0 ) ’ 


r-Q— = -\[2fc\C\zt - 


r A - + U - 1 


Z2 - 


,(0) 


d G° r 

- Cl— Zi -I-, 

or g 


/t i ( i 2) /t (0) \ 

rp( 2 ) jrp(0) 

V ) 


(28) 


(29) 


(30) 

(31) 


r 


dz 3 

dr 


= —C2 


dzi 

~dP 


- V ad V^ + V{Vad- V)zi 


+ C3 



cut) Za + 



d In l(. 0) 

—71-23 

amr 


VX7 


Z 4 + I 


(32) 


dZA 

dr 


= UV (4 -KT + qt) za-VV (k p + — )Vz 2 - VVz 3 - VV J°, 


Xp 


and equations (1B80I) and (1B81 1) as 


dz 


~^ + fC 1 B z t = l\ 1 0 /2 ^ + 


dr 


Cl 


(GYG° 

V2gci 


f r o _ dzt _ r r o _ . {G q + G q ) 
— J'-bZh —7^7 — —Y 2 fC c zi -I- —j =— — 


(33) 


(34) 

(35) 
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where the inhomogeneous terms H°, 1°, J°, Gx, and Gx with X = {r, q, q} are vectors whose k- th components are 
respectively given by , I\ k , Ji k , ^x i k , and G' x l , , which are defined in Appendix B4, and Ao and A 0 ' are diagonal 

matrices whose fc-tli diagonal components are given by A i k = Ik {Ik + 1) and A i k , respectively. Note that G^ 2 \, — G^j, = 0 
and + G^J = 0. The symbols C\, C%, C g, and denote matrices whose kj components are respectively given by 

„„J /Ak l j ,,,V, 0T . 0 1 I h 1 ] 1 , C l k l j — fd’kdj 1 . J 

u a » ( -'b > °b , ana o c , wnere o A - °o(-i)i + O oi(-i)’— °i(-i)o + L (-i)io’°c — °io(-i) T °(-i)oi’ ana tne 
definition of the coefficient CM fc is given by (IB74I) in Appendix B3. The physical quantities V , U, rA, Ti, r, /, and ao in the 
set of equations OUT) to (TT5T) are defined as 


dlnp ( °) v _ d\nM r r l _ dlnp (0) l rflnp (0) 
d In r ’ d In r ’ d In r fi dlnr 


/ <91np\ 

U^vL’ 


t = a 0 t, f 


I 47T 12 

V 3 (Jo ’ 


(36) 


and (Jo = \/ GM/R 3 , and the definition of the other physical quantities is given in Appendix Bl. 

The set of differential equations derived above are partial differential equations with r and t being the independent 
variables. Since we are interested in mean flows driven by an unstable linear oscillation mode having a complex eigenfrequency 
to, the forcing (inhomogeneous) terms in the equations are proportional to e~ 2wit , where u>i denotes the imaginary part of to. 
To make analyses simple, we look for solutions whose time dependence is also given by the factor e~ 2ullt , that is, we assume 
that the partial derivatives d/dr is replaced by the growth or damping rate 7 = —2u>i. This simplifying assumption may 
be justified when we are interested in mean flows driven by self-excited oscillation modes. Since such self-excited oscillation 
modes must be continually pumped by certain destabilizing mechanisms (e.g., opacity mechanism) even if their amplitudes 
are saturated by some mechanisms like nonlinear couplings between many different modes, we believe that we have to take 
account of the effects of this continual pumping of the oscillation modes by introducing the growth rate into the formulation. 
From equations (13411 and (13511 . we obtain 


Zh = ZnZl + Z12Z2/C1 + Z 13, 
zt = Z21Z1 + Z22Z2/C1 + Z23, 
where 

Z 11 = V2p\N^ 1 C 1 B C° c , Z 12 = —qWj^Ag/ 2 , Z 13 = 

Z 21 = V2^f\NoiC 0 c, Z 22 = /Wq/CbAJ 72 , Z 23 = /Wp/Cg, G " - 7 W 0 1 Gq + Gq 


V^gci 


Vzigci 


W 10 = / 2 CLC°b + 7 2 E, W 01 = fcWs + 7 2 E, 


and E is the unit matrix. Substituting equations (ED and EHt into equations (BUI) to E3). we finally obtain 


dz 1 
dr 

dz 


V 


_(3 + rA-^-)E + A; /2 Zn 


zi ~ 7—E - 

Xp 


V „ Ay 2 Z! 


*0 *-12 

Cl 


22 + 7 — 24 - 7 -R 0 + Aq / 2 Zi 3 + —, 
Xp cop 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 


= " c i ( AfC^r + 7E) 2 ! - 

dz 3 


rA-^+U-l + — )£ + V2fC\Z 2 2 
Ti Xp 


z 2 + —z 4 _ V2/ciCiZ 23 - + —, (43) 

Xp g 


r -7 — = —C 2 [7 (24 - V adV 2 2 ) + V ( V a d ~ V) Z l] + C 3 


dZ4 


(e T - a T ) 24 + ( e P + — ) 2 2 


dlni™ Ao , r0 

- 23 - 77W Z4 + 1 > 

d In r tv 


= -V"V ( + — ) Fz 2 - VVZ3 + VV( 4 -k t + a T ) 24 - W J u , 


Xp 

where we have used 

P (2) - V z , R o 

P (0) Xp Xp 


(44) 

(45) 

(46) 

(47) 


and the fc-th component of the vector R° is given by 

Rk = j 0 Y£Q m (p)dn, 

and the definition of the symbol Q^ 2 \p) is given by (IB34II . The set of differential equations from (1421) to (1451) is regarded as 
the mean flow equation solved in this paper. The matrix W = —jA 0 ' Z 12 is a symmetric matrix corresponding to the matrix 
W discussed, for example, by Lee & Saio (1997) in the traditional approximation. In fact, replacing 7 by kj makes W reduce 

to W. We also note that \/ 2 /C A Z 2 2 = ^A^ 2 Zh^ . 
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For the boundary conditions applied at the stellar center, we require that the functions z\ and Z 2 are regular and that 
ds/dt = 0, which leads to 


/y( 2 ) 

\^T(°) _ 


V ad V 


n( 2) 


p(°)gr 


„( 2 ) 


ro Q 


-V{V-V ad ) = - 


iu;£ • Vs' 




-7- 


Q ( 2 ) (s). 


(48) 


The outer boundary conditions applied at the stellar surface are given by the conditions dp/dt = 0 and L r = 47rr 2 ctsbT 4 , 
from which we derive 


_ p (2) u< 2) 

^ p(°)gr rao 
and 


icj£ ■ Vp' 
p(°)gr 


(49) 


r ( 2 ) 

Jj'p 


rp(2) 


where 


ST {1) = T (1) + £ r 


£ ■ VL 

r(0) 


dT (0) 

dr 


(i) 


+ 4^(^-VdT (1 ) + - 


1 S 2 T(°) 


2 dr 2 


+ 6 


fSTWY 

vw J 


+ 1^ 


,£r dT(!) 


T(°) 


and we have assumed that Pfl is constant in the outer envelope. 


(50) 


(51) 


2.3.3 Lagrangian perturbations of second-order 

Assuming that the Lagrangian displacement vector £(cc) is infinitesimal, we may expand the perturbed velocity field v = ( n ) 
at x + £ as 


(52) 


Vi(x + £) = Vi (x) + Vi- l j(x)^j(x) + -Vi;j;fc(x)£j(x)£fc(x) H- 

= vfl^x) + vflflx) + vfl\x) + v$(x)£i(x) + v$(x)&(x) + ^Vi°J ;k (x)£i(x)£k(x) + ■■■, 

where the semicolon indicates the covariant derivative, and repeated indices imply the summation from 1 to 3. The perturbed 
velocity field v(x) is expanded in terms of the oscillation amplitude, which is also assumed to be small, as 

m (x) = v- 0) (x) + vflflx) + vfl\x) H-, (53) 

where u® is the unperturbed field. Carrying out zonal averaging, we obtain to the second-order of the wave amplitude 


Vi(x + £) = vfl ] + + vW(x)£j(x) + ^v^. k {x)^(x)^ k (x), (54) 

where we have used vfl = 0. Hence, the zonally averaged Lagrangian velocity perturbation of second-order is defined as 


5v ( fl = Vi(x + £) -vf\x) = vf\x)+vfl ) j {x)^{x) + JaO£j(x)£ fc (a:), 
where 


„(!) 


j(x)£j(x) = = -Re (iivfi-jCfl . 


(55) 


(56) 


The additional terms v^flj + in equation (1551) are called Stokes corrections (or Stokes drift for velocity field). 


2.4 Angular momentum conservation in wave-mean flow interaction 

If we employ the Lagrangian mean theory of wave-meanflow interaction, in which the position vector x is divided into the 
mean x = x and the Lagrangian displacement £(cc,t), that is, x = x + £(x,t), we may obtain in the Cowling approximation 
a mean flow equation given by (see, e.g., Lee 2013; see also Grimshaw 1984) 


d ■ 


.dp' 


P S *(» + €) = - v -(€^ 

where p is the effective density, for which p = p (e.g., Brihler 2014), and the total time derivative d/dt is defined as 


d d 


d 


-1 d 


1 d 


= tt: +Vt(x + £)— + v e (x + £)-™ + «*(* + £)—- „ 
at ot or roO r sin 0 ocp 

and, if i^ 0 '* = Vg 0 ^ = 0 as assumed for uniformly rotating stars, 


(57) 


( 58 ) 
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Vr(x + £) = 8vf?\ vg(x + £) = Sv^ 2) . (59) 

In the following, we set p = p ® for simplicity (e.g., Biihler 2014). The zonally averaged specific angular momentum in the 
2 -direction at x + £ is given, to second order of perturbation amplitudes, by 


£(cc + £) = [{x + £) x v (x + £)] ■ 


(v^ + 5v^ [£e cos 6* + (r + £ r )sin6>] — 5v cos 9 — sin 9 


= *c°> + *( 2 > 


where e z is a unit vector along the 2 -axis, and 

£ (0) = (rsm9) 2 Q, 


(60) 


(61) 


£( 2 ) _ ^ t ,;„ n„,( 2 ) _j_ ^ A,/ 1 ) < 


: r sin 9v^ + r sin + 


sin 9 + - C</> v 0 1> ) cos # + [(£r sin 9 + & cos 9 ) 2 + 12, 


and 


= ni 15 + v ( °]£j for a = r, 9, (p. 

The mean flow equation m is now given, correct to second order of perturbation amplitudes, by 


d 

dt' 


^ + [ Sv &>° + 

or 


d 5v„ 2 ' > d 


r 89 


e m = 


,( 0 ) 


v- € 


. dp 1 

9<p J 


(62) 


(63) 


(64) 


Making use of equations (16111 and (1621) . we rewrite equation (16411 as 


r sin 9- 


dv 


T) 

<t>;j 


1 p' dp' 


& -Cm ..crri ai 3“ 


dt '* J ' p(°) p(°) dip 


dv i J) 


dt 




dv 


( 1 ) 


dt 


sin 9 + 


dV ^\e-^ 9V ^ 


dt 


dt 


cos 9 


+212 


(£ r sin 9 + £e cos 9) sin 9 + Vg 1 ^ cos 9^ + + 2rf2 £j sin 2 9 + Vg£j sin 9 cos 9 j 




(65) 


(2 s ) (2 s ) (2') 

where we have used, to eliminate the second order perturbations vi- , Vg , and , the ^-component of the equation of 
motion: 


sin6 §i v f ) + 2nsin6 ( cos 9vf ) + sin0u< 2) ) = -sinflt^ 1 ’ + 


( 66 ) 


The left-hand-side of (1651) is given by the sum of products of first order perturbations associated with the oscillation mode. 
Using the perturbed continuity equation p 1 = —V • ( p ®£) and the ^-component of the equation of motion for first order 
perturbations, we rewrite the term (p^) -2 p'dp'/dip in equation (l65l) as 


1 p' dp' 


_rm _tni oj. _(o) ^ ’ ( £ oj. ) + £ ’ ^ 1 nj. 


p(°) pi°l dip 


.dp' 


1 dip 


1 dp' 


p(°l d(p 


= 


. dp' 
’ dip 


-«-v 


r sin 9 -M - 2 r sin 9Q ^sin 9v < J' > + cos 9vg^ ^ 


(67) 


which we substitute into equation (1651) to prove the identity. 

We may use the meanflow equation (1651) or (1641) along with equation (1661) to check numerical consistency. To simplify the 
computation, we integrate equation (IM1) over spherical surface to obtain 


d_ 

dt 


^l’ ( ' 2 ' > ^ + 2rf2 (jjvfp 2 ^ sin 2 9 + Svg 2 ’ 1 sin 9 cos 9^ ^ 
where 

W{r) =irr 2 Iia((Cp')) , 
and 


d 


27rr 2 p(°) dr 


W(r), 


n 7T n2 

( f) = / 
Jo Jo 


f sin OdOdcj). 


( 68 ) 


(69) 


(70) 
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Figure 1. dW for the prograde l = \m\ = 2 p-mode of the /3 Cephei star model (solid line) and for the retrograde V = m = 2 r 36 -mode 
(dashed line) and the prograde l = \m\ = 2 p 3 o-mode (dotted line) of the SPB star model, where we assumed £2 = 0.1 for the p- and 
p-modes and Cl = 0.4 for the r-mode. We plot 5 x dW for the p-mode. 


The function W(r) may be regarded as a work function (e.g., Unno et al 1989), and dW/dr > 0 and dW/dr < 0 respectively 
indicate the excitation and damping regions for the oscillation mode. We use equation m, in stead of (RHl) , to see the 
numerical consistency. 

For later use, we rewrite (1681) into a non-dimensional form C = 77, where 


C = cilPy + 2cifl 



ra o 


sin# + 


SjP 

rao 




(71) 


n=™(rd- lm CjL- 

2 \ or \ r pgr 


+ U + 2 ) Im (T-P- 
alnr / \ r pgr 


= mdW, 


(72) 


and 


(4 2) ) = <V 2) ) /r 2 0.. (73) 

It may be useful to plot the function dW, as defined by equation dZSD, to indicate the locations of mode excitation and 
damping regions in the outer envelope of the stars for the oscillation modes discussed in this paper. Figure 1 shows dW for the 
prograde l = \m\ = 2 g 3 o-mode and retrograde l' = m = 2 r 36 -mode of a 6 Mq SPB star model and for a prograde l = \m\ = 2 
p-mode of a 15 Mq /3 Cephei model, where we assume S2 = 0.1 for the g- and p-modes and 0. = 0.4 for the r-mode, and the 
physical parameters of the SPB star and /3 Cephei star models are given in §3.1 and §3.2, respectively. As shown by the figure, 
the strong excitation occurs in the region at r/R ~ 0.95 and a damping takes place at r/R ~ 0.92 below the excitation region. 


3 APPLICATIONS TO PULSATING MASSIVE STARS 

To show an application of the above formalism to rotating massive stars, we apply it to rotating, slowly pulsating B (SPB) 
stars and /3 Cephei stars. As background equilibrium models for oscillation calculation, we use stellar models computed with a 
standard stellar evolution code, where no effects of rotational deformation are considered in evolution calculation. The opacity 
tables used for evolution and oscillation calculations are those computed by Iglesias & Rogers (1996). In this paper, we are 
interested in axisymmetric mean flows driven by non-axisymmetric oscillation modes of rotating stars, where the oscillation 
modes are excited by the K-mechanism associated with the iron opacity bump. Non-adiabatic oscillations of uniformly rotating 
stars are computed with the method employed in Lee & Saio (1987). The quantity / in that latter paper is set to 1 so that 
the effects of the centrifugal acceleration and of the corresponding rotational deformation of the equilibrium structure are not 
taken into account in our computations of the oscillations (see e.g., Lee & Baraffe 1995). 
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Figure 2. Top: second-order perturbations us.k = ^] k /R a 0t u H,k = and ux : k = \ /R&0 ;KS a function of r/R for the 

even retrograde g§ 4 -rnode of / = |m| = 2 of a 6 Mq main sequence star model at Q = 0 . 1 , where the solid lines, dash-dotted lines, dashed 
lines, long dashed lines, and dotted lines represent respectively the components corresponding to k = 1, 3, 5, 7, and 9. Bottom: density 
( 2 ) ( 2 ) ( 2 ) 

maps of the functions u). (r,6) / Rctq, (r,Q)/Ri to, and v\ '(r, 9 )/Rctq for the mode. The functions are normalized by their maximum 

amplitudes. Red and blue colors correspond to positive and negative values, respectively. 


3.1 Slowly Pulsating B-type Stars 

In SPB stars, low frequency g- and r-modes are excited by the opacity bump mechanism (Dziembowski et al. 1993; Gautschy 
& Saio 1993; Lee 2006; Aprilia et al. 2011). Here, we calculate mean flows driven by unstable low frequency g- and r-modes 
of a 6 Mq main-sequence star, whose physical parameters are typical of a SPB star: log (L/Lq) = 3.2328, logT e ff = 4.1982, 
R/Rq = 5.55, and X c = 0.1237, X — 0.7, and Z = 0.02. Because good numerical consistency in the sense discussed in §2.4 
is attained only for Q <0.1 for the g-modes, we restrict our discussion of g-modes to the case of B = 0.1. Since most of the 
unstable p-modes of the SPB model have frequencies |w| B for S2 ~ 0.1, we can use the expansion length j max ~ 5 to obtain 
sufficiently accurate eigenfunctions, and we set fc ma x = 10 for second order perturbation calculations. 


3.1.1 g-modes 

Let us first discuss mean flows driven by low frequency g-modes. As numerically shown by Aprilia et al. (2011), both prograde 
and retrograde low frequency g-modes are excited by the opacity bump mechanism in rapidly rotating SPB stars. Most of 
l = \m\ retrograde g-modes of the star, that would be unstable if the star were non-rotating, are stabilized by rapid rotation 
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as a result of their coupling with stable l > \m\ p-modes. As a consequence, only a small number of them survive as unstable 
modes in rapidly rotating SPB stars. This is not the case of prograde l = \m\ g-modes and most of them remain unstable 
even at rapid rotation. 

In Figure 2, we plot, in the top panels, the expansion coefficients us,k = flg\ k /R(J q, UH,k = Vfl\ k /R<Jo, and ur,k = 
vlg\ k /Rao as a function of r/R for the even retrograde ( 730 -mode of l = \m\ = 2 for S2 = 0.1, where the normalized eigenfre- 
quency uj = w/cro is id = (0.4672, —5.7 x 1CP 6 ). Here, the amplitude normalization of the linear mode is given by Si 1 = 1 at 
the surface. Note that l = \m\ implies the mode is an even mode. As shown by Figure 1, the functions us,k, UH,k, and Ur,k 
have large amplitudes only in the outer envelope layers. Note also that the peaks of the functions at r/R ~ 0.95 correspond 
to the place at which the mode excitation strongly occurs. The reason for this behavior is that, in addition to the fact that 
mode amplitudes can be large near the surface because of low density, thermal diffusion grows below the surface of upper 
main-sequence stars. Therefore, deposition/extraction of momentum may take place in these layers (e.g., Lee et al. 2014). 
These wave-mean flow interactions thus drive both zonal and meridional flows there (see e.g. Mathis et al. 2013 in the case of 
low-mass stars). The amplitude of the zonal component Ur,k is much larger than those of the meridional ones, us,k and UH,k, 
by several orders of magnitude. Therefore, the kinetic energy of the induced azimuthal rotational motion is larger than that 
of generated meridional circulation. 

The reason that the azimuthal flow is much larger than the radial/latitudinal flows can be simply understood from the 
propagation of the mode pattern. A purely adiabatic non-axisymmetric mode is a standing wave in both the radial and 
latitudinal directions. However, the mode propagates around the equator of the star and can be thought of as a traveling wave 
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in the azimuthal direction. Consequently, the mean flows created by the wave are much larger in the azimuthal direction, 
because the traveling wave propagates in one direction (either prograde or retrograde). In the radial/latitudinal direction, 
the standing wave, which is composed of traveling waves propagating in opposite directions, contains components that cancel 
each other out such that there is no net flow in these directions. In other words, since the modes contain angular momentum 
only in the z-direction, the flows they induce are only in the azimuthal direction. For a weakly non-adiabatic mode such as 
those discussed here, the modes also produce non-zero mean flows in the radial/latitudinal directions. However, the azimuthal 
flows will still be much larger. 

Figure 2 also shows for the 530 -mode the density maps of the components v ^ (r, 8)/Rao, Vg 2 \r, 8)/Ra o, and ( r , 8)/Rao, 
in the bottom panels, where these functions, which depend on both radial distance and latitude, are normalized by their 
maximum amplitudes. The vertical axis is the rotation axis (x = y = 0), and the equatorial plane is given by z = 0. The 
radial and the azimuthal velocity components (t)f* (r, 8)/Rao and v^\r, 8)/Rao) are symmetric about the equator, while the 
latitudinal one (vg 2 \r,8) / Rao) is antisymmetric. As shown by the maps for the retrograde <730 mode, the pulsation driven 
mean flows near the surface have the velocity components Vg ; toward the equatorial plane and the component v) ' heading 
inwards. The <f> component 1 M 1 is positive, indicating acceleration of the rotation near the equatorial surface. As suggested 
by the density maps, the amplitudes of pulsation driven mean flows tend to be confined in the equatorial region. This reflects 
the properties of low frequency modes of rotating stars, that is, their amplitudes are also confined in the equatorial region, 
particularly for rapidly rotating case (see, e.g., Berthomieu et al 1978; Bildsten et al 1996, Lee & Saio 1997; Townsend 2003ab). 

We see that pulsations induce a differential rotation 5Q(r,8 ) = (r, 9)/(r sin#), which is a function of both radial 

distance and latitude. This is the signature of pulsation-driven transport of angular momentum both in the vertical and in 

the horizontal directions (e.g., Mathis 2009). It is the first time that such wave-driven differential rotation is computed in 2-D. 

Indeed, previous computations treated reduced cases of the radial shellular rotation (e.g., Talon & Charbonnel 2005; Mathis 

et al. 2013) or focused on the dynamics in the equatorial plane (e.g., Rogers et al. 2013). Like in previous computations 

performed for low-mass stars by Talon & Charbonnel (2005) and Mathis et al. (2013), we see in massive stars that the 

meridional circulation can be multicellular both in radius and in latitude. These patterns differ from properties of meridional 

circulations driven by large-scale differential rotation in models where waves are not taken into account. Such a computation 

is of interest for development of 2-D models of rotating stars (Rieutord 2006), in which waves must be taken into account in a 

near future. In the case of the 530 -mode, the amplitudes of the functions are confined in the equatorial regions. In the context 

(2) 

of active massive stars, such as Be stars, it is interesting to note that for the retrograde mode, the velocity perturbations v^ 
at the surface becomes positive in the narrow equatorial region, while v) ; is negative. The positive zonal component can then 
help the surface layers to reach the critical angular velocity needed to eject matter in the circumstellar environment. 

In Figure 3, we show, in the top panels, the expansion coefficients us,k , UH,k, and Ur,k for the even prograde 530 -mode of 
l = \m\ = 2 at Cl = 0.1, for which Co = (0.4025, —7.9 x 10 -6 ). The corresponding 2-D density maps of the velocity components 
Vr 2 \r,8)/Rao, Vg 2 \r, 9) /Rao, and v^\r,9)/ Rao are given in the bottom panels in the same figure. As in the previous case, 
the functions us,k, UH.k, and Ur,k have large amplitudes only in the outer surface layers. As in the previously studied case of 
the retrograde 530 -mode, the amplitude of Ux,k is much larger than those of us,k and UH,k- From these density maps, the flow 
patterns driven by the prograde 530 -mode in the surface equatorial region is opposite to those by the retrograde 530 mode, 
that is, the radial velocity v) 1 is towards the surface and the components Vg ; show flows out of the equatorial plane. This 
different behavior between retrograde and prograde 5 -modes shows that it is necessary to sum over all the modes that are 
excited in a given star to be able to understand the total wave-driven zonal and meridional mean flows and their net effect. 
This is the reason why r- and p- modes have also to be considered. 


3.1.2 r-modes 

r-modes are rotationally induced retrograde modes and constitute a subclass of inertial modes, for which the Coriolis acceler¬ 
ation is the restoring force. In the slow rotation limit (i.e., fl —>■ 0), r-modes associated with the degree l' and the azimuthal 
index m have an asymptotic co-rotating frame frequency given by u = 2 mfl/ [l'(l' + 1 )], and the displacement vector £ is 
dominated by its toroidal component (i Tp). Since numerous l' = \m\ r-modes of odd parity are also excited by the opacity 
bump mechanism in SPB stars (e.g., Lee 2006), it is necessary to examine how r-modes drive mean flows in rotating SPB 
stars. 

The top panels of Figure 4 show the expansion coefficients us,k, UH,k, and ux,k as a function of r/R, while the bottom 
panels of the figure give the corresponding 2-D density maps of the three components vi 2 ~ > (r, 8)/Rao, Vg 2 \r,9)/Rao, and 
( r , 8)/Rao for the odd r 36 -mode of l' = m = 2 at fl = 0.4, for which Co = (0.2199, —2.1 x 10 -6 ). As in the case of g-modes, 
the amplitudes of the expansion coefficients become large in the outer envelope. Moreover, the zonal flow (ut, k) is still much 
larger than the meridional one (given by us,k and UH,k). As shown by the 2-D density maps, gross flow patterns driven by 
the r-mode are similar to those by the retrograde 530 -mode, although the zonal acceleration near the surface is much more 
prominent than that for the retrograde 5 -mode. 
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Figure 4. Same as Figure 2 but for the l' = \m\ = 2 odd r 36 -mode at fl = 0.4. 


3.1.3 Lagrange velocity perturbations 

The vertical Lagrangian velocity perturbation 5v can be an important quantity when we consider radial transport of gaseous 

matter in the interior of stars. In Figure 5, we show the 2-D density maps of 5vfl /Roo for the even retrograde and prograde 

P 30 -modes of l — \m\ = 2 at Cl = 0.1 and for the odd r 36 -mode of l' = \m\ = 2 at Q — 0.4. Comparing the bottom left panels 

( 2 ) 

of Figures 2, 3, and 4 with the corresponding ones in Figure 5, we see that there appear no essential differences between vi- 
(2) 

and Svi ' for the low frequency modes. When zonal flows in the surface equatorial layers are accelerated by the retrograde 
modes, the radial flows are driven inwards. 


3.2 /? Cephei stars 

In addition to exciting g- and r-modes in cool B stars (SPB stars), the iron opacity bump mechanism also excites p-modes 
in hotter B stars (/3 Cephei stars). We calculate mean flows driven by unstable low radial order p-modes of a 15Mq main- 
sequence model, whose physical parameters are typical of a /3 Cep star: log (L/L®) = 4.5928, logT e g = 4.394, R/R® = 10.8, 
and X c = 0.130. For this model, l = |m| low radial oder modes with normalized frequency u) between 2 and 5 become 
pulsationally unstable. Since the method of calculation used in this paper for p-modes of rotating stars is not necessarily 
appropriate for very rapid rotation Cl ~ 1, as suggested by full 2-D computations of p-modes in rapidly rotating stars 
(Lignieres et al. 2006; Reese et al. 2006), we restrict the present discussion to the case of slow rotation rates Cl ~ 0.1. 
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Figure 5. Density maps of the functions 5v^ ( r , Q)/rcro , from left to right panels, for the retrograde ( 730 -mode and the prograde ( 730 -mode 
of l = \m\ = 2 at Cl = 0.1, and for the l' = \m\ = 2 r 36 -mode at Q = 0.4 for the 6 Mq main sequence model. The functions are normalized 
by their maximum amplitudes. Red and blue colors correspond to positive and negative values respectively. 


Figure 6 and 7 are respectively for the retrograde l = m = 2 p-mode of w = (2.584, —7.9 x 10 -7 ) and for the prograde 
l = \m\ = 2 p-mode of to = (2.634, —4.5 x 10 -7 ). Note that, if we count the number of radial nodes of the eigenfunction Sq 
in the way described in Unno et al (1989), we obtain k = n p — n g = —3, suggesting that the mode should be classified as 
a <? 3 -mode, where n g and n p are the numbers of p-type and p-type nodes of the eigenfunction S q . However, we simply call 
these modes a p-mode because it behaves as a p-mode in the envelope. 

As shown by Figures 6 and 7, and as in the case of SPB stars studied above, the zonal component ( Ur,k ) is dominant over 
meridional ones (us,k and un,k)- The flow patterns v^ and v^ in the surface layers are quite similar between the retrograde 
and prograde p-modes, which is contrary to the case of low frequency p-modes. This is because the Coriolis terms proportional 
to are not important to determine these velocity components. Zonal acceleration in the surface equatorial region occurs 

for the retrograde mode, although deceleration takes place for the prograde mode. The Lagrangian velocity perturbations 
Svi- 2 ^ of the p-modes are depicted in Figure 8, which shows that the additional terms can be significant enough, that is, 

even if vi 2 ^ is positive in the surface equatorial region, Svi 2 ^ becomes negative. 


4 DISCUSSION AND CONCLUSION 

As a consistency cheque of our mean flow computations, in Figure 9 we plot the functions £ and 77 for the prograde p- and 
p-modes, and the retrograde r-modes depicted in Figure 1. Note that equation ([Mji is well satisfied for these modes. The 
two functions agree very well for the p- and p-modes, but they show disagreement in the outer most layers for the r 36 -mode. 
Although the reason of the disagreement for the r-mode is not well understood, possible reasons might be that the number 
of radial nodes of the eigenfunctions is large, which makes it difficult to correctly compute the eigenfunctions, and that the 
toroidal components of the displacement vector are significantly dominant over the the other components although equation 
(16611 is determined by the phase differences between the less dominant eigenfunctions. 

Assuming an initial uniform rotation, we have derived the governing equations for zonal and meridional axisymmetric 
mean flows driven by unstable non-axisymmetric oscillation modes in rapidly rotating massive main-sequence stars, where the 
magnitude of mean flows are assumed to be of second-order of the oscillation amplitudes. The governing equations are a set of 
coupled linear ordinary differential equations for the second-order quantities with inhomogeneous terms, coming from products 
of the eigenfunctions of linear oscillation modes. To demonstrate the applicability of the formalism and its importance for 
astrophysical studies, we have computed zonal and meridional axisymmetric mean flows driven by non-axisymmetric p- and 
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Figure 6. Same as Figure 2 but for the l = 


\m\ = 2 even retrograde p-mode of the 15Mg main sequence star model where Q = 0.1. 


r-modes in SPB stars and p-modes in j3 Cephei stars, where oscillation modes are assumed to be excited by the iron opacity 
bump mechanism. 

For most of the oscillation modes considered in this paper, these mean flows have large amplitudes only in the surface 

regions of the stars, particularly in the regions where the oscillation modes are excited. The first interesting point to note 

(2) 

is that for low frequency retrograde g-modes and r-modes excited by the opacity bump mechanism, v ^ ' can be positive in 
the surface equatorial regions, indicating that surface fluid elements could be accelerated in the same direction as the surface 
stellar rotation. The velocity fields generated in the surface layers are related to the Reynolds stresses of the waves and 
transported fluxes that drive exchanges between waves (oscillations) and zonal and meridional mean flows (e.g., Bretherton 
1969; Andrews & McIntyre 1978a; Holton 1982; Mathis et al. 2013; Belkacem et al. 2015a). For the zonal component, we 
provide for the first time the 2-D geometry of the wave-driven differential rotation given by v^\r, 6)/ (r sin 0ao). This provides 
us with information on the transport of angular momentum by waves both in the vertical and in the latitudinal directions 
(e.g., Andrews & McIntyre 1978a; Mathis 2009). 

Moreover, we computed for the first time 2-D wave-driven meridional circulation. It can have multi-cellular pattern in 
the radial direction, while its latitude-dependence depends on the studied modes. If we use Eq. (1) to roughly estimate the 
magnitude of the rotationally-driven meridional circulation, we have ti r; MC/ R&o ~ 10 -7 x fl 2 for a 10 Mq ZAMS star. This 
number suggests that pulsation-driven mean meridional flows have comparable magnitude to rotationally-driven meridional 
circulation if the oscillation modes have amplitudes Sq ~ 10~ 4 to 10~ 3 at the surface, depending on the mode. To estimate 
the magnitudes of velocity fields of pulsation-driven mean flows, we need to theoretically estimate the amplitudes of pulsation, 
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Figure 7. Same as Figure 6 but for the l = |m| = 2 even prograde p-mode. 


which has always been a difficult problem. One possible way to estimate oscillation amplitudes for the k —mechanism is to 
employ a weak non-linear theory of oscillations (see e.g., Lee 2012), although it is another very difficult task to apply this 
theory to oscillation modes in rapidly rotating stars. As suggested by Lee (2012), the amplitudes of this magnitude Si 1 ~ 10” 4 
to ICR 3 is probably attainable in SPB stars. This is the reason why wave-driven mean zonal and meridional mean flows must 
be taken into account when studying the evolution of rotating massive stars, as it is done for low-mass stars (e.g., Talon & 
Charbonnel 2005; Mathis et al. 2013). 

In this paper, we have applied several simplifying assumptions to make the formulation and calculation of mean flows 
tractable. First, we used the radiative transfer equation to derive the governing equation for the second-order variables, even 
for convective regions in the interior. This crude treatment may be justified for massive main-sequence stars having only 
weak surface convection layers, but may not be justified for low-mass stars with thick surface convection zones. Note that 
K-driven pulsations do not have any contribution to mean flow in the convective core of massive stars, unless the effect of the 
viscous force on inertial modes is taken into account. However, the induced transport of momentum is negligible compared to 
transport by convective eddies (Browning et al. 2004). In addition, we ignored the effects of centrifugal force on the equilibrium 
structure, oscillation calculation, and mean flow computation for rapidly rotating stars. This neglect of the centrifugal force 
effects may not be a serious drawback for the mean flows driven by low frequency g- and r-modes. However, these effects have 
to be taken into account for mean flows driven by p-modes in stars rotating as rapidly as Cl ~ 1 (not treated in this paper). 
Next, we assumed an initial uniform rotation, and neglected the possible existence of critical layers and breaking regions, 
at which waves suffer strong dissipation and efficient exchange of momentum between mean flows and waves may take place 
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Figure 9. C (dashed line) and 7Z (solid line) for the 930 -mode and r 36 -mode of the SPB star model and for the p-mode of the 0 Cephei 
star model, from left to right panels, where we use H = 0.1 for the g- and p-modes and Q = 0.4 for the r-mode. 


(Rogers et al. 2013; Alvan et al. 2013). These processes will be studied in the near future. Finally, low-frequency modes could 
also be excited stochastically by turbulent convection in the core of massive stars (e.g., Rogers et al. 2013; Lee et al. 2014; 
Mathis et al. 2014). In particular, stochastically excited waves appear to be very important for the ejections of matter by Be 
stars (Neiner et al. 2012; Lee et al. 2014). The theory developed in this work will be applied to these stochastically-excited 
modes in the future. 

This work shows that pulsation-driven flows can be as significant as rotation-driven flows in pulsating massive stars. 
Therefore, it is important to take them into account when studying the rotational and chemical evolution of massive stars. In 
addition, pulsation-driven flows can transport angular momentum to the surface layers, resulting in acceleration or deceleration 
of rotation velocity of the fluid of stars. This angular momentum transport mechanism might help to form a circumstellar 
gaseous disc around Be stars, although the magnitudes and directions of the mean flows in the radial direction suggested in 
this paper are not necessarily favorable for the mechanism to be viable. In this paper, we have treated the oscillations as 
standing waves in the radial direction. If we assume low frequency modes become progressive in the surface layers because of 
their very low frequencies as suggested by Ishimatsu & Shibahashi (2013), the way of transport of angular momentum by the 
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waves could be different from that by standing waves. This possibility will be pursued in a future paper concerning angular 
momentum transport by stochastically excited oscillations of massive stars. 
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APPENDIX A: SPIN-WEIGHTED SPHERICAL HARMONIC FUNCTIONS 

Spin-weighted spherical harmonic functions s Y; m (#,<)>) may be defined as eigenfunctions of the differential equation 

as.yr = -(i+s)(i-s+i) 3 Yr, (ai) 


or 


m s Y t m = -(l-s)(l + s + l) s Yr, 

where the differential operators 9 and 9 are defined for s Y ; m with —l^s^las (Newman & Penrose 1966) 


(A2) 



(A3) 



(A4) 


For a given value of s, the function s Y; m is normalized as 
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J (sY™ 1 )* a Y™ 2 sin 6 Wdcf> = 5i 1 i 2 5m 1 m 2 , 
where 

(.V7T = (-i) s+m - s ’S'r m - 

We note that oY™ = Y™, and it is convenient to use 

So yt = -VXi -1 y, m , 9-1 yt = VMoYr, ss 0 yr = -MoYT 
9 oYT = VAliYT, 9S 0 Yr = -AioYT, 

3 -iFi° = V20Y?, SiW 0 - 0 , S-iYi° = 0 , SiY® - -V20Y?, 


(A5) 

(A 6 ) 

(A7) 

(A 8 ) 

(A9) 


where A* =1(1 + 1). 

Angular integration of a product of three spin-weighted spherical harmonic functions can be evaluated by using the 
formula given by 


■ h 

h 

h 


mi 

m2 

m 3 

- / 

. s i 

S2 

S3 . 



/ v m l 


Y™ 2 V' m 3 JQ 
; 2 2 i 2 S3U-, 


= 47T 


47T 


2/1 + 1 / 2/2 + 1 / 2/3 + 1 7 Zi 


47T 


47T 


^2 


I3 \ f h h 1 3 


— Sl —S2 —S2 / V m l m 2 m 2 


(A10) 


^2 ^3 


where dfl = sin 9d6dd>, and 

\m 1 m2 m2, 

zero values only when mi + m 2 + m 3 = 0 , si + S2 + S 3 = 0 , and h + I2 ^ I3 ^ \h — Zi|. 


is Wigner 3 -j symbol (Rose 1957; Edmonds 1968). This integration has none 


APPENDIX B: SECOND-ORDER EQUATIONS FOR MEAN FLOWS AND SPIN-WEIGHTED 
SPHERICAL HARMONICS 

B1 Introducing the basis associated to Spin-Weighted Spherical Harmonics 

We introduce a new set of basis vectors e q and Gq defined by 
_ G.6 + i Gq, _ gq — ie^, 

for which 

Gj. * Gq — Gf • Gq — 0 ; G q • Gq — Gq • Gq — 0 , Gq • Gq — 1 . 

Using e r , G q and Gq, we may rewrite the displacement vector £(x) as 
£(*) = &e r + £ q G q + iqGq, 


(Bl) 


(B2) 

(B3) 


where 

4? ’ 45 s/2 ' 

The differential operator V may be rewritten as 
„ d d d 

V — G r 7- b Gq -j— + Gq -j— , 

or dq dq 

where 
d 

n — 
dq 

d 

f)- — 

dq 

We find that 

<9 v m _ 

o-oR ~ ~ 

oq 

^ _ 

— — 

dq 


(B4) 


(B5) 


= J_( 

s/2r V 

d . 1 d\ 

dd 1 sin 9 d(j>) ’ 

(B 6 ) 

1 / 

d 1 d \ 

(B7) 

_ V2r 1 

dd 1 sin 8 dij>) 


m, 1 rr~ \+rn 

= - 7=-V Mill , 

V 2 r 

(B 8 ) 


m 1 /~T~ -\7-rn 

= Wr^~ lh ■ 

(B9) 
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Using the basis vectors e q and e q and spin-weighted spherical harmonics sY™, we rearrange the series expansions m 
to (11311 as 


0 = 1 

Jmax 


Jmax 

fr = E Ot 3 °^7 - iT^Sol^) 


3=1 


and similarly the expansions (O to (EH) as 


'“max 

,,(2) - E 410*1°, 

k =1 

1 felIlaX , V 

f E (41^° +4 2 ^ou ; ?) 


'“max 

4 2) = -^ E (4V<< -4 2 1'3ou ; ?) . 


(BIO) 

(BH) 

(B12) 

(B13) 

(B14) 

(B15) 


If we write v^ = irf-iY°e q +irfiY°e q using / = ^/ 47 t /3 9 and sin# = ^Sit/SiYi = — -\/ 87 r/ 3 -iUi 0 , the Coriolis term 
in equation (1191) reduces to 

v (2) • Vu (0) + u (0) • Vu (2) = 2i/ [- (w^iY ! 0 + uf-ilf) e r + (u^-iYi 0 - u^oH 0 ) e, + (v^iY ) 0 + rf’o^ 0 ) ej , (B16) 

and hence the r, q, and g-components of the momentum equation m are written as 


dv 


( 2 ) 


dt 


1 dpV pW g(2) 

i - y „(o) v ~ T u ’ 


an, 


( 2 ) 


at 


+ 21 / - .for. 0 ) + = G ?’. 




( 2 ) 


at 

where 


+ 2i/ (n^i^ 0 + 4 2) oU°) + 


dq 

1 ^P (2) _ r ( 2 ) 
p<°> ag ? ’ 


G ^ = —v' • Vt>' + g 


P' 

0 ( 0 ) 


e r + -^y-J^-Vp' = G^e r + G^fl e q + G q ^e q . 


We rewrite eauations (1201). (IB17I). (1221). and the radial comDonent of eauation m into a non-dimensional form: 


^ a Vr 2S> _ a p *- 2 - 1 ^ v ^ 

dr ra o p(°) rcro 


- (3 + 


dlnp (0) \ n^ 2) tf (2) 


dlnr y roo p(°)cro’ 


a p® 


= — + 2 i Cl / ( ] - 


( 2 ) 


ar p(°)pr 3r rcro 


,( 2 ) 


w 0 , % 


( 2 ) 


rcro 


rcro 


d In p^gr p ^ p*- 2 ^ Gr 2 '* 


dlnr p(°)gr p(°) g ’ 


(B17) 

(B18) 

(B19) 

(B20) 

(B21) 

(B22) 


( 2 ) 


^ a l: 

ar £,(°) 


a t (2) 


-C 2 


+C3 


a r (2) „ a p (2) ” (2) 




3r T(°) “ dr p(°) rcro 


c,, + 


1 \ p 1 


( 2 ) 

^o7 


Xp/ P 

rp( 2) 


+ (er — ar) 


rp(2) 

T( °) 


dlni™L$. 2) 1 


+ / (2) , 


dlnr p/o) i/V " T(°) 


1 \ p' 


= ^V(4- h T + ar) ^ - UV ( k, + - ) ^ - UV^ - UVJ 


( 2 ) 


,( 2 ) 


/( 2 ) 


Xp/ P' 


o(0) 


ar r(°) 

where 

cro = y/ GM/R 3 , r = a 0 t, Cl = fl/a 0 , f = f/a 0 , 

Ll 0) = 4nr 2 Ffl\ L (2) = Anr 2 Ffl\ F r (0) = -A ( 0 ) dT ( 0 ) /dr, 


do) 


(B23) 

(B24) 

(B25) 

(B26) 
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Cl = 


°0 


v = - 


(r/B) 3 

g/r Mr/M ’ 
dlnp^ 


C 2 = 


4tt r 3 c p 


d In r ’ 
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( d In k \ 

KT - \d^f); 
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(B27) 
, (B28) 

p 

(B29) 

(B30) 


J (2) = c 3 


p e 
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p(°) e(°) CpUQ 


v' • Vs' + 


i+(JL + -C\ + 
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+ fi 

Si 
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+” ; —J 
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vv 


J (2) = Q (2) (k) + Q (2) (p ) - C (2) + ( 
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1 d T' 
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T(°) J 
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2 h dp 2 \p(°) 


2 hdT 2 \T(°)J 


<o) T (o) 1 dh, p' T' 

1 h dpdT p(o) T(°) ] 


and h denotes functions that depend on p and T, vff = v q 2 ^e q + v q 2) e q , and 


Vh = r e, 


5 q dq +e? dqJ ' 


(B31) 

(B32) 

(B33) 

(B34) 

(B35) 


Note that equations (IB18I) and (IB19I) will be used to derive algebraic equations relating the variables v q 2 ^ and v q 2 ^ to the 
variables and p^. 


B2 Expansion of non-linear terms in terms of Spin-Weighted Spherical Harmonics 

In this Appendix, we give explicit expressions for various products of linear wave functions, which are given by series expansion 
in terms of spin-weighted spherical harmonic functions s Y/ m in the basis coordinates ( r,q,q ) and basis unit vectors e r , e q 


To evaluate the terms v' ■ V»', we need covariant derivatives of the velocity vector, or the displacement vector ( see 
Schenk et al 2002) 
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(B43) 
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(B44) 


& = E F -™-iYi 

l^\m\ 


where we have used the notation for the covariant derivatives VkP , and 



Jmax 

fo m (r) = rJ2 S h (r)S l i., 
3=1 



and 


= A; (A, - 2). 

Thus, the terms v' ■ V«' in equation (IB20I) may be given as 


(B45) 

(B46) 

(B47) 

(B48) 

(B49) 

(B50) 

(B51) 

(B52) 

(B53) 

(B54) 

(B55) 

(B56) 

(B57) 



Since the term p'Vp' in equation (IB20I) is given by 
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7W=\ EMHR ^ (o^ m o^) ( 0 y-™_ lF -) e , - ( 0 y-- +1 y-) e? |, (B59) 

5 M 2 ^(-i) m K [(4)* df‘> {oY^oY,™) - (/L\) V\ ( + 1 y, 7 m _ 1 y J ™) - ( 4 \) (_ 1 y l - m + 1 y l ™)] 


we obtain 
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-^Ec- 1 )"* 
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The term V • ( p'v') in equation (IB30II may be evaluated as 


v-(pv) = b^E 


1 9 


Ph ( r 2 Q r r M 2 Hl 2 ) +‘S' i 2 r Ar ,Pi 
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The terms VjjA' • V hT' and £ • Vs' in equation (IB31I) may be given by 

\ I* rj-if 
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We also note that the Coriolis term and the term Vp® in equation (1191) are given as 
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B3 Projection of second-order equations onto Spin-Weighted Spherical Harmonics sY™ 

Multiplying equations (IB21I) to (IB24I) by oY k , and carrying out angular integration over spherical surface, we obtain 
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(B70) 
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The symbol [• ■ •] in equation (IB74I) has been defined in Appendix A. 

Multiplying equations (IB 181) and (IB19I) by i Y k and -iY k , and carrying out angular integration over spherical surface, 
we obtain 
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We may obtain a set of algebraic equations using equations (IB75I) and (IB76I) . The sum (IB75I) + (IB76I) gives 
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where we have used 
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B4 Explicit spectral coefficients of non-linear terms in the basis of Spin-Weighted Spherical Harmonics 

Here we give explicit expressions for the quantities G®, G^ k , H^\ lj?\ J k 2 \ which are projections onto spin-weighted 
spherical harmonic functions: 
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To evaluate Q we need to calculate second derivatives of thermal quantities. For example, for the specific entropy s 
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For the density p, we have 
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For the opacity, which is given as a function of p and T, we have 
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